# TODO: Add comment
# 
# Author: guochun
###############################################################################

source("./sourceAll.R")
#aggregated species
data=read.table("./objects/data/bci5.txt",header=T)
sp1=data[data$sp==unique(data$sp)[3],]
sp2=data[data$sp==unique(data$sp)[220],]
par(mfrow=c(1,2),mai=c(0.9,0.9,0.3,0.2))
plot(x=sp1$gx,y=sp1$gy,xlab="x",ylab="y",cex=0.8,col=rgb(62/255,93/255,119/255))
plot(x=sp2$gx,y=sp2$gy,xlab="x",ylab="",cex=0.8,col=rgb(62/255,93/255,119/255))
unique(data$sp)[3]
unique(data$sp)[220]

#pattern changes
library(spatstat)
X <- rThomas(100, 0.02, 10)

plot(X,main="")
Z <- as.im(function(x,y){ 5 * exp(2 * x^2 - 1) }, owin())
Z$v=Z$v/max(Z$v)
Y=rthin(X,Z)
plot(Z,ribbon=FALSE,main="")
plot(Y,main="",add=T)

#real habitat heterogeneity
ph=getData("PH.dbf")
ph=im(ph)
plot(ph,main="")
ph$v=(ph$v/max(ph$v))^15
X2=rpoispp(0.15,win=c(ph$xrange,ph$yrange))
Y2=rthin(X2,ph)
plot(Y2,add=T)
X <- rThomas(0.0006, 10, 40,win=owin(ph$xrange,ph$yrange))
plot(X,main="")


ph=getData("PH.dbf")
ph=im(ph)
#windows()
plot(ph,main="")
ph$v=(ph$v/max(ph$v))^5
X3=rthin(X,ph)
plot(X3,add=TRUE)

k=im(getData("K.dbf"))
mg=im(getData("Mg.dbf"))
ph=im(getData("PH.dbf"))
plot(k,main="",ribbon=FALSE)
plot(mg,main="",ribbon=FALSE)
plot(ph,main="",ribbon=FALSE)


